#pragma once
#include "elisa/elisa.h"
#include <iostream>
#include <sstream>
#include <algorithm>
#include <string>
#include <vector>
#include <gaol/gaol>

#include "Modeler.h"

#ifdef ELISA_PERFORMER
#include "elisa_chrono.h"
#endif

using namespace std;
using namespace elisa;

/* Problem: filter design problem

@unpublished{BiniMourrain1996,
author    = {D. Bini and B. Mourrain},
title     = {{Handbook} of {Polynomial} {Systems}},
month     = {November},
year      = {1996}
}

128 real roots in the box [-1.0e8,1.0e8]^9
*/
class Filter1 : public Modeler
{
private:
	Model _model;

	RealVar a;
	RealVar b;
	RealVar m[7];
	RealConstraint c1;
	RealConstraint c2;
	RealConstraint c3;
	RealConstraint c4;
	RealConstraint c5;
	RealConstraint c6;
	RealConstraint c7;
	RealConstraint c8;
	RealConstraint c9;

public:
	Filter1() // constraints might be wrong..
		: a(Interval(-1.0e8,1.0e8), "a"),
		  b(Interval(-1.0e8,1.0e8), "b"),
		  c1(0.01 == m[1]*m[3]*m[5]),
		  c2(7/500 == a*m[3]*b),
		  c3(2/25 == pow(a,2) + pow(m[0],2)),
		  c4(37/50 == pow(b,2) + pow(m[6],2)),
		  c5(0.9401 == pow(m[1],2) + pow(m[2],2) + pow(m[3],2) + pow(m[4],2) + pow(m[5],2)),
		  c6(0.038589 == pow(m[1],2)*(pow(m[3],2) + pow(m[4],2) + pow(m[5],2)) +
		     pow(m[2],2)*(pow(m[4],2) + pow(m[5],2)) + pow(m[3],2)*pow(m[5],2)),
		  c7(-0.0081 == (m[4]*m[6] - m[5]*b)*(m[0]*m[2] - a*m[1])),
		  c8(m[3]*(m[0]*m[1]*m[2]*b + a*m[4]*m[5]*m[6]) - (7/500)*(pow(m[5],2)+pow(m[1],2)) == 39/25000),
		  c9(2.7173/4 == pow(m[6],2)*((9401/10000) - pow(m[5],2)) +
		     pow(b,2)*((9401/10000) - pow(m[3],2) - pow(m[4],2)) - 2*b*m[4]*m[5]*m[6]),
		  _model()
	{
		for(int i=0; i<7; i++) {
			string name("m["), e("]");
			ostringstream oss;
			oss << (i+1);
			name += oss.str();
			name += e;
			new( m + i ) RealVar(Interval(-1,1), name);
		}
	}

	elisa::Model & model()
	{
		_model = c1 && c2 && c3 && c4 && c5 && c6 && c7 && c8 && c9;
		return _model;
	}

	void print()
	{
		cout << '\t' << a.getName() << " = " << a.domain() << endl;
		cout << '\t' << b.getName() << " = " << b.domain() << endl;
		for(int i=0; i<7; i++) {
			cout << '\t' << m[i].getName() << " = " << m[i].domain() << endl;
		}
	}

	RealVar & varX()
	{
		return a;
	}
	RealVar & varY()
	{
		return b;
	}
};
